#include "fermiqcd.h"

//
// Action and metropolis are implemented here:
// (please check signs and formulas)
//

void generate_next_config(mdp_field<float> &phi, 
			  coefficients& param, 
			  int nhits=1) {
  // unpack parameters from param object
  float sigma=sqrt(param["alpha_1"]);
  float alpha_2=param["alpha_2"];
  float alpha_3=param["alpha_3"];
  float alpha_4=param["alpha_4"];
  // get a reference to lattice object
  mdp_lattice &lattice=phi.lattice();
  // auxiliary variables
  mdp_site k(lattice);
  float x, y, c, deltaSprime;

  // first do even sites, then do odd sites, k
  for(int parity=EVEN; parity<=ODD; parity++) {    
    forallsitesofparity(k,parity) {
      // on each site try nhits new variables
      for(int hit=0; hit<nhits; hit++) {
	// generate x with required gaussian distribution
	x=lattice.random(k).gaussian()*sigma;
	// y is old value of phi(k)
	y=phi(k);
	// compute variation in action if phi(k)=y -> phi'(k)=x
	c=phi(k+0)+phi(k-0)+
	  phi(k+1)+phi(k-1)+
	  phi(k+2)+phi(k-2)+
	  phi(k+3)+phi(k-3)+
	  phi(k+4)+phi(k-4);
	// the term that multiplies alpha_3 is P (CHECK SIGN)
	deltaSprime=	  
	  -c*(x-y)+
	  alpha_2*(pow(x,4)-pow(y,4))+
	  alpha_3*((x*x-alpha_4)/pow(x*x+alpha_4,2)-
		   (y*y-alpha_4)/pow(y*y+alpha_4,2));
	// cout << x << " " << y << " " << c << " " << deltaSprime << endl;
	
	// metropolis accept/reject here
	// cout << "deltaSprime=" << deltaSprime << endl;
	if(exp(-deltaSprime)>lattice.random(k).plain()) {
	  phi(k)=x;
	  break; // if accept stop multi-hit and go to other site
	}
      }	      
    }
    // perform communication if in parallel
    phi.update(parity); 
  }
}

float average(mdp_field<float>& phi) {
  float avg=0;
  mdp_site k(phi.lattice());
  forallsites(k) avg=avg+fabs(phi(k));
  mdp.add(avg);
  return avg/phi.lattice().global_volume();    
}

float compute_one_over_F(mdp_field<float>& phi, 
			     float alpha_3, 
			     float alpha_4) {
  float P, one_over_F=0;
  mdp_site k(phi.lattice());
  forallsites(k) {
    P=alpha_3*(phi(k)*phi(k)-alpha_4)/pow(phi(k)*phi(k)+alpha_4,2);
    one_over_F+=exp(-P);			      
  }
  mdp.add(one_over_F);
  return one_over_F/phi.lattice().global_volume();
}

//
// This program starts from a cold configurations phi(k)=0
// and generates nconfig sets (save as dataset000001, dataset000002, ...)
// on each dataset phi it computes and prints
// O1=(sum_k phi(x))^2 and O2=(sum_k phi(k))^4
//
// The program assumes average of O1 and O2 over nconfig(s) is done off-line
//

int main(int argv, char** argc) {
  // open communication channels if in parallel
  mdp.open_wormholes(argv,argc);
  coefficients param;
  // ////////////////////////////////////////////////////////
  // PHYSICAL PARAMETERS OF COMPUTATION
  // I am making them up now
  // You probably want to input them from file or keyboard
  // ////////////////////////////////////////////////////////
  // lattice ndim 
  int ndim=(int) val(prompt("", "LATTICE_NDIM", "5"));
  // lattice size  
  int n=(int) val(prompt("", "LATTICE_SIZE", "6"));
  // number of configs to be genetared
  int nconfig=(int) val(prompt("", "NUMBER_CONFIGS", "100"));
  // number of configs to be skipped
  int nthermalize=(int) val(prompt("", "NUMBER_THERMALIZE", "10"));
  // number of hits in multihit metropolis
  int nhits=1;            
  // 5D lattice
  int box[5]={n,n,n,n,n}; 

  // IMPORTANT: this is the input for A(a)/(2 a^2)
  param["F"]=1;           // this is one because is has to be determined
  // input alpha_1 = 20+a*m(a)^2
  param["alpha_1"]=val(prompt("", "ALPHA_1", "30")); 
  // input alpha_2 = lambda(a)/a
  param["alpha_2"]=val(prompt("", "ALPHA_2", "10")); 
  // this is zero because it has to be determined
  param["alpha_3"]=0.0; 
  // input alpha_3 = A(a)/(2*a^2)  
  float alpha_3=val(prompt("", "ALPHA_3", "0.1")); 
  // input alpha_4 = B(a)*a^3
  param["alpha_4"]=val(prompt("", "ALPHA_4", "1.0")); 
  // ////////////////////////////////////////////////////////

  // declare lattice object
  mdp_lattice lattice(ndim,box,default_partitioning0,torus_topology);
  // declare field phi(k)
  mdp_field<float> phi(lattice); // main field
  // declare site object k of coordinates k(0), k(1), ..., k(4)
  mdp_site k(lattice); 

  // auxiliary variables
  char filename[128];
  JackBoot jb(nconfig,1);
  
  // START FROM COLD CONFIGURATION
  forallsites(k) phi(k)=0;
  // perform communication if in parallel
  phi.update(); 

  // do a number of steps to thermalize the metropolis
  for(int config=0; config<nthermalize; config++) {
    generate_next_config(phi,param,nhits);
    cout << "mean_k |phi(k)|=" << average(phi) << endl;
  }

  // loop over all field configurations
  for(int config=0; config<nconfig; config++) {
    // call metropolis and compute next
    generate_next_config(phi,param,nhits);  
    // generate filename and sava dataset for reuse
    sprintf(filename,"dataset_tune_%.6i",config);
    phi.save(filename);
    cout << "mean_k |phi(k)|=" << average(phi) << endl;

    // compute observable <exp(-P)>
    jb(config,0)=compute_one_over_F(phi,alpha_3,param["alpha_4"]);

    // print observable for each config
    cout << config << "\texp(-P)=" << jb(config,0) << "\n";
  }

  // print mean and bootstrap error over configurations
  jb.plain(0);
  cout << "1/F = <exp(-P)) = " << jb.mean() << " +/- " << jb.b_err() << "\n";

  // close communication channels
  mdp.close_wormholes(); 
  return 0;
}







